# Political Surveys Bias Voters' Self-Reported Economic Perceptions
# Script 8: National Economic Perceptions by EU Vote Results Figure

# Jack Bailey
# University of Manchester
# jack.bailey@manchester.ac.uk


# 1. Housekeeping ---------------------------------------------------------

# The following code requires that you have first opened the associated
# project file. If not, click the button on the top right and open it
# ("Open Project...").

# Note also that the Bayesian methods I use below rely on having installed
# the probabilistic programming language Stan. If you have yet to install
# Stan, please see https://mc-stan.org/users/interfaces/


# Set random seed

set.seed(666)


# Load packages

library(tidyverse)
library(tidybayes)
library(jbmisc)     # https://github.com/jackobailey/jbmisc
library(brms)
library(here)


# Load model

m3 <- readRDS(here("_output", "m3.rds"))



# 2. Triptych plot --------------------------------------------------------

# We're going to create three plots that describe how the treatment affected
# the national retrospective economic perceptions that our subjects reported
# based on who they voted for at the 2017 referendum. We're also going to
# respect the ordinal nature of the data too. First, we'll create a plot that
# shows how the responses of incumbents in the treatment group compared to
# those in the control group.

# Use m3 to predict 'e' where 'r' = 'Leave' and 't' = 'Treatment'

inc_t <- 
  add_fitted_draws(m3,
                   newdata = tibble(r = "Leave", t = "Treatment")) %>% 
  ungroup()


# Use m3 to predict 'e' where 'r' = 'Leave' and 't' = 'Control'

inc_c <-
  add_fitted_draws(m3,
                   newdata = tibble(r = "Leave", t = "Control")) %>% 
  ungroup()


# Compute inc_t - inc_c to determine the average treatment effect for each
# response option conditional on having voted for the incumbent. Then we'll
# pick and rename the variables we'll need to plot it.

inc_cate <-
  inc_t %>%
  mutate(cate = .value - inc_c$.value) %>% 
  select(r, t, resp = .category, cate)


# Use m3 to predict 'e' where 'r' = 'Remain' and 't' = 'Treatment'

opp_t <- 
  add_fitted_draws(m3,
                   newdata = tibble(r = "Remain", t = "Treatment")) %>% 
  ungroup()


# Use m3 to predict 'e' where 'r' = 'Remain' and 't' = 'Control'

opp_c <-
  add_fitted_draws(m3,
                   newdata = tibble(r = "Remain", t = "Control")) %>% 
  ungroup()


# Compute opp_t - opp_c to determine the average treatment effect for each
# response option conditional on having voted for the incumbent. Then we'll
# pick and rename the variables we'll need to plot it.

opp_cate <-
  opp_t %>%
  mutate(cate = .value - opp_c$.value) %>% 
  select(r, t, resp = .category, cate)


# Again, we do the same as above, but for those who did not vote at the 2017
# referendum. These instructions are the same as above, but substitute "Non"
# for "Opp".

# Use m3 to predict 'e' where 'r' = 'Nonvoter' and 't' = 'Treatment'

non_t <- 
  add_fitted_draws(m3,
                   newdata = tibble(r = "Nonvoter", t = "Treatment")) %>% 
  ungroup()


# Use m3 to predict 'e' where 'r' = 'Nonvoter' and 't' = 'Control'

non_c <-
  add_fitted_draws(m3,
                   newdata = tibble(r = "Nonvoter", t = "Control")) %>% 
  ungroup()


# Compute non_t - non_c to determine the average treatment effect for each
# response option conditional on having voted for the incumbent. Then we'll
# pick and rename the variables we'll need to plot it.

non_cate <-
  non_t %>%
  mutate(cate = .value - non_c$.value) %>% 
  select(r, t, resp = .category, cate)


# Now that we've computed the estimated CATEs for those respondents who voted
# for the incumbent, an opposition party, or did not vote at the last election,
# we need to combine the estimates into a single dataset. We'll also convert
# the 'p' column to a factor at the same time so that we can control the order
# in which its options appear in our plots.

cates <- 
  rbind(inc_cate, opp_cate, non_cate) %>% 
  mutate(
    r =
      r %>% 
      factor(levels = c("Leave", "Remain", "Nonvoter"),
             labels = 
               c("Leave",
                 "Remain",
                 "Nonvoter")),
    resp =
      resp %>% 
      factor(labels =
               c("Lot Worse",
                 "Little Worse",
                 "Stayed Same",
                 "Little Better",
                 "Lot Better")
      )
  )


# Now we'll save the CATEs to disk so that we can include them in the paper.

saveRDS(cates, here("_output", "cates_m3.rds"))


# We'll also use these data to create summaries that we can use as labels

cate_labs <-
  cates %>% 
  group_by(r, resp) %>% 
  summarise(
    median = median(cate),
    .groups = "drop")


# Now, we have the data for the plot and the labels too, so can begin to plot
# the data. First, we tell ggplot2 that we want to use the 'cates' dataset we
# just created, that we want it to plot each option in the 'p' variable in a
# different facet, and that we want to add a dotted line at the zero mark.

eu_fig <- 
  cates %>% 
  ggplot(aes(x = cate, y = resp)) +
  facet_wrap(~ r) +
  geom_vline(xintercept = 0,
             linetype = "12",
             colour = bailey_colours("grey6"),
             size = .5)


# Next, we'll add our data to the plots and also the labels that we just created
# too. At the same time, we'll also add in points that show the median value of
# each distribution.

eu_fig <- 
  eu_fig +
  stat_halfeye(aes(slab_fill = r,
                   slab_colour = r),
               .width = .95,
               point_size = 1,
               interval_size = 0.5,
               slab_alpha = 0.7,
               slab_size = 0.5,
               height = .7,
               normalize = "xy",
               show.legend = F) +
  geom_point(data = cate_labs,
             aes(x = median,
                 y = resp), colour = "white",
             size = 0.3, show.legend = F) +
  geom_text(data = cate_labs,
            aes(label = scales::percent(median,
                                        accuracy = 0.1,
                                        suffix = ""),
                x = median + .002,
                y = resp),
            family = "Cabin",
            size = 2,
            vjust = 2,
            hjust = 0.7,
            inherit.aes = FALSE)


# Finally, we need to make the plot look nicer.

eu_fig <- 
  eu_fig +
  scale_y_discrete(labels = function(x) str_wrap(x, width = 8)) +
  scale_x_continuous(labels = scales::percent_format(accuracy = 1, suffix = ""),
                     breaks = seq(-.15, .15, by = .05)) +
  scale_colour_manual(values = bailey_colours("blue", "red", "grey"),
                      aesthetics = "slab_colour") +
  scale_fill_manual(values = bailey_colours("blue8", "red8", "grey8"),
                    aesthetics = "slab_fill") +
  coord_cartesian(xlim = c(-.15, .15)) +
  theme_bailey() +
  theme(axis.title.y = element_blank(),
        axis.ticks.x = element_line(lineend = "round"),
        axis.text.y = element_text(hjust = 0.5)) +
  labs(x = "Conditional Average Treatment Effect (Percentage Points)")



# 3. Replication details --------------------------------------------------

# Save session information

save_info(here("_output", "_session_info", "009_eu_fig.txt"))


# One last thing...

thanks()
